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We generalize the screened Korringa-Kohn-Rostoker (SKKR) method for solving the correspond¬ 
ing Kohn-Sham-Bogoliubov-de Gennes (KSBdG) equations for surfaces and interfaces. As an ap¬ 
plication of the newly developed theory we study the quasiparticle spectrum of Au overlayers on 
a Nb(100) host. We find that, within the superconducting gap region, the quasiparticle spectrum 
consists of Andreev bound states (ABS) with a dispersion which is closely connected to the underly¬ 
ing electronic structure of the over layer. We also find that the spectrum has a strongly k-dependent 
induced gap. The properties of the gap is discussed in relation to the thickness of the overlayer, and 
it is shown that certain states do not participate in the Andreev scattering process. 


I. INTRODUCTION 

The theory of Bardeen, Cooper, and Schrieffer (BCS) 
successfully describes the universal properties of conven¬ 
tional (s-wave) superconductors 1 , but it can not be ap¬ 
plied easily to inhomogeneous systems where the wave 
number, k , is not a good quantum number. The gen¬ 
eralization of the well-known Hartee-Fock method with 
the introduction of the concept of mixed particle-hole ex¬ 
cited stated 2 3 -l yields the Bogoliubov-de Gennes (BdG) 
equation^. In this description the standard momen¬ 
tum operators are replaced by field operators, which 
have the advantage that they are able to describe in¬ 
homogeneous systems. However, this is only a mean- 
field theory and can not be considered as a predictive 
approach to allow the computation of material-specific 
properties. For that purpose a density functional the¬ 
ory (DFT) was constructed by Oliveira, Gross and Kohn 
(OGK) 5 . In this theory the ground state energy is proved 
to be a unique functional of the p(f) charge density and 
the x(r) = (Ati(r)) anomalous density. Later, the 
concept is further developed into a multicomponent den¬ 
sity functional theory for the combined system of elec¬ 
trons and nuclei (phonons) 6 8 . The usefulness of the 
OGK approaclPhas been demonstrated by Suvasini et al. 
where they introduced a simple semi-phenomenological 
parametrization of the exchange-correlation functional 9 . 
Despite the simplicity of this approximation, they were 
able to describe many features of the bulk niobium in 
the superconducting state, which are accessible to exper¬ 
iments. Using this semi-phenomenological parametriza¬ 
tion, one can derive a set of equations which allows the 
self-consistent solution of the following coupled Kohn- 
Sham-Bogoliubov-de Gennes (KSBdG) eigenvalue equa¬ 
tions in atomic (Rydberg) units 9 


(H e (r) ~ Ef) u n {r) + A eff (r)v n (r) = e n u n (r ), (la) 
(H e (r) - E f ) v n (r) - A * eff (r)u n (f) = -e n v n (r), (lb) 


where H e (r) = — V 2 + U e //(r) is the single-particle 
Hamiltonian, and the wavefunction is decomposed into 
an electron-like part u n (r) and a hole-like part v n (r). 
The effective electrostatic and pairing potentials are 



A eff(r) = A x(r), (2b) 

where V ext (r) is the external potential (e.g. the Coulomb 
attraction from the protons). The p{f) charge and y(r) 
anomalous densities can be calculated from the wavefunc¬ 
tion components: 

p(r) = 2 X [\ u n(.r)\ 2 f(tn) + K(r)| 2 (l - /(£«))] (3a) 

n 

X(r) = Ty„(r)<(f)(l - 2/(e„)). (3b) 

n 

f(e n ) is the Fermi-Dirac distribution function, E® c [n\ is 
the usual exchange correlation energy for normal elec¬ 
trons and A is a semi-phenomenological adjustable pa¬ 
rameter (it can be site-dependent). It should be noted 
that the zero point of the energy scale is the Fermi level. 

The past few years have shown a growing interest in the 
study of superconductor based heterostructures 10 12 . It 
is known that such inhomogeneities in the pairing poten¬ 
tial can lead to bound quasiparticle states. These states 
have been found theoretically in superconductor - nor¬ 
mal metal - su percon ductor heterostructures 13 ^ and also 
in other systems 14 15 . The Andreev reflection 16 1 has been 
identified as the key effect which results in such bound 
states, called Andreev bound states (ABS): an electron, 
with energy lying in the superconducting gap, arriving 
from the normal metal to the superconductor - normal 
metal (S/N) interface is retro-reflected as a hole and a 
Cooper pair is formed in the superconductor. While a 
great many theoretical works were dedicated to study 
the Andreev reflection and the ABSPtHl it was done 
on model systems only, their material specific dispersion, 
their ’’band structure” has never been calculated (nor 
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observed experimentally) to date. Within the framework 
of a tight-binding model, Annett and coworkers also in¬ 
vestigated such heterostructures 22 23 l They have shown 
the existence of ABS within the gap, and pointed out 
effects associated with the interplay of the gap and the 
normal van Hove peaks 22 . The next logical step in this 
series of investigations is first principles calculations for 
real materials. In this paper we address this problem by 
developing a multiple scattering theory (MST) for the so¬ 
lution of the KSBdG Eqs. 0 for surfaces and interfaces 
of S/N heterostructures. As presented in Refs. 9]and-17. 
the application of constant pairing potentials gives a very 
good estimation of the self-consistent solutions. This fa¬ 
cilitates to model a S/N system with a finite, constant 
pairing potential on the superconducting host and zero 
on the normal metal over layers. In order to treat this 
semi-infinite geometry, a Green function based method 
is needed, like for examp le a Screened Korringa-Kohn- 
Rostoker (SKKR) met ho within MST. 

The present paper is organized as follows. In Sec. II 
we generalize the SKKR method for the solution of the 
KSBdG Eqs. 0- Sec. Ill is denoted to the computa¬ 
tional details. In Sec. IV we illustrate the power of the 
developed method by studying the quasiparticle ’’band” 
structure of the niobium (Nb) - gold (Au) system. Fi¬ 
nally Sec. V is devoted to the summary. Some technical 
details are provided in the Appendix. 


II. GENERALIZATION OF THE MST FOR 
SUPERCONDUCTORS: THE BDG-SKKR 
METHOD 


The central problem of the DFT calculations is the 
solution of the KSBdG Eqs. 0 in order to determine 
the single-particle wavefunctions and the corresponding 
eigenvalues. However, the single-particle Green function 
contains all information about the ground state. The 
local DOS, the anomalous and charge densities can be 
obtained from the single-particle Green function. Conse¬ 
quently, if the single-particle Green function is obtained, 
it is not necessary to calculate the Kohn-Sham orbitals. 
In this section we show how the SKKR method can 
be generalized to get the single-particle Green function 
for multilayered systems corresponding to the KSBdG 
Eqs. 0 . Here we do not try to follow every single step 
of the derivation, as it would be too extensive for this pa¬ 
per, instead we just show how the most important quan¬ 
tities, concepts and formulas needs to be modified due 
to the presence of holes. Most interim derivations can 
be performed in analogy of the normal state MST, well 
described for example in Ref. [26[ The first step in this 
generalization is to decompose the BdG Hamiltonian in 
the following way: 


( 4 ) 


where 


Hoif) 


/_V 2 -E f 0 \ 

V o v 2 + e f )' 


( 5 ) 


V(r) 


WffW -Vefftf) ■ 


(6) 


In the KKR method, the potential is usually treated 
in the so called muffin-tin approximation, ie. the po¬ 
tential is written as a sum of single-domain potentials 
centered around each lattice site, n, namely V e // (r) = 
Vn{r). It is usually assumed, that V n {r) is spher¬ 
ically symmetric. This is not a necessary assumption 
for the theory, however, MST for a general shape of po¬ 
tentials are still not common even for the normal state. 
Therefore, in what follows we still restrict ourselves to 
spherical atomic potentials. In our approach, we as¬ 
sume the same form for the effective pair interaction as 
well, namely A e ff(r) = J2 n The potentials V n {r) 

are of muffin-tin type which means that V n {r) = 0 and 
A n (t ) = 0 if r = \r n \ > *S n , where S n is the muffin-tin 
radius. 


A. Operator formalism and the free-particle Green 
function 

The resolvent of the BdG Hamiltonian can be defined 
as 


g(z) = (zI-H BdG )-\ (7) 

which has the usual property 

Q{z*) = g oy. (8) 

At the real axis the up- and down-side limits of Q(z) are 
defined by 

G(z = e±iO) — G ± (e). (9) 

Similarly to the normal state, we can define the T- 
operator as 

T(z)=V + VG(z)V = J2 Tnm ( z )> ( 40 ) 

n,m 


where r nrn (z) is the so called scattering path operator 
(SPO) which comprises all possible scattering events be¬ 
tween the sites n and m, including now the Andreev re¬ 
flection as well. Since V is Hermitean, 

T{z*) = T{z)K ( 11 ) 

The two different, generalized eigenfunctions of HsdG 
can be obtained from the Lippmann-Schwinger equation: 

= <p(e) + ^(e)r ± (e)v?(e), 


Head?) = H 0 (r) + V(r), 


( 12 ) 
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where ip(e) is a generalized, multicomponent eigenfunc¬ 
tion of Hq and Qo(z) is the resolvent corresponding to Hq. 
Here we emphasize the difference from the normal state 
now, that is the wavefunctions also have a hole part. 

Following the footsteps of normal state MST, we define 
the following orthogonal and complete basis set: 

^ ( “ ri j , (13a) 


’ < 13b > 

where L = (/,ra) is the usual composite index, and 
jl( z ,r) = il{p e r)Y L (r), jt(z,r) = ji(p h r)Y L (f), (14) 
ji(x ) is the spherical Bessel function of the first type, 

p e = yjE F + z, p h = \]E f - z. (15) 


In a first step the free-particle Green function is derived 
which is related to the structure constant describing the 
structural properties of the investigated system. Using 
the definition of the basis abobe, Eqs. (13), in terms of 


contour integrations - commonly used in MST - it can 
be shown, that the Green function of free particles has 
the following form: 

Go b (z,r,r')= 6 ab H l{ z X>)\jl{ z X<)] X (16) 


where 

H a L (z,r n ) = -i p a h a L (z,f n ) = -ip a h?(p a r n )Y L (r n ). 

(17) 

Here we used the notation r< = min (r, r'), r> = 
max (r, r f ) with the choice of Im{p e } > 0 and Imjp^} > 
0, a, b denote the electron, hole indices, and hf(x) = 
( x ) is the Hankel function of the first type, while 
h%(x) = — h+(x). We defined the x operator for any 
arbitrary function, /, as fl(z,r n ) x = fi(p a r n )Y L (r n y . 
Therefore, it is clear, that Gg 6 is diagonal in indeces a 
and 6, and the hole part of the free-particle Green func¬ 
tion can be obtained from its electronic part 

G h 0 h (z; f, r 1 ) = -G e 0 e (-z; r, r 1 ). (18) 


structure calculations, a scalar relativistic generalization 
of the BdG theory is needed. To arrive to such a theory, 
one needs to start from the relativistic Dirac-Bogoliubov- 
de Gennes (DBdG) equations already worked out in the 
literature 2 ^^. An analogous scalar relativistic form of 
the BdG equations can be obtained quite straightfor¬ 
wardly by neglecting not only the spin-orbit coupling 
term but all relativistic corrections to the pairing po¬ 
tential as well. By suppressing the explicit dependence 
on the complex energy, on a log-scale (x = logr) these 
coupled equations for the radial part - since both the po¬ 
tential and the pairing potential is spherically symmetric 
- can be written as follows: 


-^Qi ( x ) = ~Qt( x ) + Uf(x)P?(x) + e x A (x)P l h (x), 

(19a) 

= P?(x) + e x B e (x)Qf(x), (19b) 

Or) = -Q?{x) + Upx)P h (x) - e x A*(x)P?(x), 

(19c) 

^-Ppx) = Ppx) + e x B h (x)Q?(x), (19d) 


where the wavefunctions are defined as 

p t( x )\ =p x (M x )\ 
Pi h (x)J \Vl(x)J 


and 


C,‘W = |btil+e“(V'W+*), 

B-(x) = 1 + izZM, 

c z 

B h (x) = 1 - Z + 


( 20 ) 

(21a) 

(21b) 

(21c) 

(21d) 


In MST, the scattering matrices and the scattering solu¬ 
tions are obtained by matching the solutions of the above 
equations inside the muffin-tin sphere to the solutions 
outside. 


C. Single-site scattering 


B. Scalar relativistic Bogoliubov-de Gennes 
equations 

Nowadays almost all electronic structure codes are 
built around what is called the ” scalar relativistic” ap¬ 
proximation, where the mass-velocity and Darwin terms 
are properly taken into consideration, but the spin-orbit 
coupling is neglected. Consequently, to be able to thor¬ 
oughly compare our results with normal state electronic 


After performing the necessary integrations in the 
Lippmann-Schwinger equation (12) and using the basis 


defined in Eqs. (13) together with the definition of the 


free-particle Green function Eq. (16), two different scat¬ 


tering solutions can be obtained and written in the fol¬ 
lowing matrix-form 


R n L ’ ab (z, r n ) = j a L (z, ? n )S ab +J2 H b A z , r n )tppz), (22) 
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where is the single-site t-matrix which is diagonal 

in L, 1/ indices for potentials with spherical symmetry. 
Equation ( [22| implies that an electron (hole) like solution 
to the Lippmann-Schwinger equation may have a hole 
(electron) like component as well. If the incoming wave 
is electron like, the solution outside the muffin-tin sphere 
can be written as 


sent at ion, 

{t n £${z)S nm }, (26) 

{G n Q m L ^(z)(l-Snm)}, (27) 

{T™’ ab (z)}, (28) 


t(z) = 
Go (2) = 

r(z) = 


R e ( ^ _ (ji(p e r) - i p e tf e (z)h+(p e r)\ . 

K L( Z ' r >-{ ip h t^(z)h+(p h r ) 

(23a) 

and if the incoming wave is hole like, it can be expressed 


R h L (z,r) 



YlmiO,</>). 


(23b) 

It should be noted that these equations are different not 
only in the electron-hole components, but also in the ap¬ 
propriate energy dependence as well through p e and p h . 
As mentioned earlier, the t-matrix can be obtained by 
matching the outside scattering solutions and the regu¬ 
lar solutions inside the muffin-tin sphere at the boundary, 
which is described in more details in the Appendix. 

Using the particle-hole symmetry®, it can be easily 
proved that 


t? e (-z) = tf\z), (24a) 

tr(-z) = -t1\z). (24b) 

These symmetry relations are independent from the 

actual form of A (r) and V(r) in the superconducting 

muffin-tin sphere. 


D. Multi-site scattering, generalized 
Faulkner-Stocks formula 


where the t(z) SPO can be determined from the single¬ 
site t-matrix and the real space structure constant simi¬ 
larly to normal state in the supermatrix formalism 

T ( z ) = (t(z) -1 — G 0 (x)) _1 . (29) 

For periodic systems in KKR theory it is possible to 
write the above equations in reciprocal space, which al¬ 
lows to obtain the SPO as a function of the wave num¬ 
ber and the energy. Finding the poles of the SPO as a 
function of k and e gives the electronic band structure. 
Butler described a one-dimensional version of KKRP^, 
which is often used as a testbed for new ideas within the 
theory. This has been done in Ref. [30| where a one¬ 
dimensional model of the Bogoliubov-de Gennes-KKR 
theory has been presented. However, since translational 
invariance is broken at the interface, to be able to cal¬ 
culate physical properties on surfaces and interfaces, we 
follow the derivation of a full real-space Green function, 
and make use of two-dimensional periodicity later. 

In the normal state MST, it had been shown by 
Faulkner and StockJ 3 ^ that the Green function can be 
obtained from the scattering path operator and from the 
scattering solutions. The derivation can be followed step 
by step, and the result for the site-diagonal part of the 
Green function can be written in an analogous form as 
well. First, the full Green function was evaluated for the 
case of | v n | > S n , | r m | > *SVn: 


A rather convenient property of the KKR method that 
it allows a transparent decoupling of the potential (de¬ 
scribed by its scattering matrix) and the structural prop¬ 
erties of the system of scattering centers (atoms). Simi¬ 
larly to the normal casd^, using the well-known expan¬ 
sion of plane waves into spherical Bessel function and 
spherical harmonics (Bauer’s identity 26 ) and the defini¬ 
tion of the free-particle Green function Eq. ( p~6| ), the free, 
real space structure constants to describe the structural 
properties can be constructed as: 


G(z,r,r) = Z(z,r)r(z)Z(z,r) x - Z(z, r)J(z, r) x , (30) 

where we used the following matrix notation (F = 
Z, or J): 

F (z\ r) = {F n,a6 (z-r)} 

= fr b (z;ri, fr b (z-,ri, •••]}, 

(31a) 


G o™L b ( z ) s a,b4n^2i L L ' L " Hl„(z,R nm )CL L ", 

L" 

(25) 

where R nrn is the vector pointing from site n to site m 
and Cf' L „ are the usual Gaunt-coefficients 26 . 

The scattering matrices, the matrices of the structure 
constant and the scattering path operator can be intro¬ 
duced in a quasiparticle-site-angular momentum repre- 


and the corresponding adjoint vector, 


F (z; r) x = {F"’ ab ( 2 ;r) x } 


/ 

\fi’ ab ( z ;r) x 1 



fr b 


< 

fs’ ab ( z ', r) x 


< 




(31b) 
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where 


J? b (z,r) = j?(z,r)5 ab , (32a) 

Z?\z, r) = ]T Ri C {z,r) [t^] C \z), (32b) 

C 

Z?\z,r) = [tT'rWRf&r), (32c) 

C 

Rf(z, r) = J? b (z, r ) + H?{z, r)tf(z), (32d) 

R?\z, r) = J? b (z, r) + tf b (z)H[(z, r). (32e) 


To calculate physical quantities, we have to continue 
the Green function inside the muffin-tin spheres by us¬ 
ing the solutions of the scalar relativistic BdG Eqs. (19) 

I 


Piif) = ~~ [°° de/(e) [ 3 k 
71 J -OO J BZ 

Xi(r) -T J de(l - 2/(e)) 


III. COMPUTATIONAL DETAILS 


In this section, we describe the technical details of the 
calculation of quasiparticle spectrum for a real supercon¬ 
ducting heterostructure using the BdG-SKKR method 
outlined in Sec. II. 

The geometry of our system builds up from two- 
dimensional translational invariant layers. The system 
comprises three regions: (i) semi-infinite bulk (Nb); (ii) 
the interface region that - in our case - consists of six 
superconducting layers (Nb), various number of normal 
metal layers (Au) and three layers of empty spheres; (iii) 
and semi-infinite vacuum. The Nb has the body-centered 
cubic (BCC) crystal structure with a lattice parameter 
a = 3.3 A. Here we do not try to investigate the effect of 
matching different lattice structures on the quasiparticle 
spectrum. Thus, for simplicity we assume BCC epitaxial 
growth for the Au overlayers and the Nb/Au BCC(100) 
heterostructure will be investigated. 

As we mentioned in the introduction, we do not 
calculate self-consistently the A e ff(r) pairing-potential, 
only the normal state calculation is performed self- 
consistently to obtain the V e ff(r) effective potential. We 
do this to simplify our first calculations, and because it 
has been shown in Ref. 0 that a good guess of the self- 
consistent pairing-potential can be the following average 



where Vws is the volume of Wigner-Seitz cell. Conse¬ 
quently, we treat the A averaged pairing potential as an 


as described in details in the Appendix. The formulas 
given above can be applied to surfaces and interfaces 
quite straightforwardly following the idea of the so called 
Screened KKR (SKKR) formalism described in Refs. l24l 
and ESI In this formalism, a special reference system is 
used to obtain structure constants that are localized in 
real space. In the supermatrix formalism we used above, 
the screening transformation can be written in a way that 
is formally exactly the same as it was presented in Ref.l25l 
Thus the whole formalism can be derived for layered sys¬ 
tems with two-dimensional periodicity and applied as the 
SKKR method prescribes. 

To perform fully self-consistent calculations for S/N 
systems, it is necessary to calculate the charge density 
and the anomalous density for layer /, which can be ob¬ 
tained from the Green function: 


ImTr a , L G a L b lV + (e,f,k I,,), 


(33a) 



Im Tr L G h L e ^’ + (e,r,k „). 


(33b) 


adjustable parameter. Since A is the experimentally ob¬ 
served gap (the gap is measured from the Fermi level), 
in principle, it should be set to equal the experimen¬ 
tal value 32 . However, with this value of the A (orders 
of magnitude smaller (meV) than the electronic energies 
(eV) involved in a normal state band structure calcula¬ 
tion), many layers are necessary to see its effect on the 
bands crossing (not just near the Fermi level), which sig¬ 
nificantly increases the computational time. Therefore, 
a model A is used here to explore the quasiparticle spec¬ 
trum. The conclusions we draw, however, does not de¬ 
pend on the size of the A parameter. 

Similarly to normal state electronic structure calcula¬ 
tions, single-site t-matrices are obtained for each layer, 
where the A averaged pair interactions can be different 
on each layer, just as the atomic potentials. In our model 
a finite, constant A pairing potential is assumed on the 
Nb layers and A = 0 Ry on the Au over layers. 


In practice, we obtain the t-matrix and the wavefunc- 
tions in the following way: first, the radial scalar rel¬ 
ativistic BdG Eqs. (19) are integrated outwards up to 
the radius of the muffin-tin sphere with different start¬ 
ing values to obtain the R^ 6 (e, r). The matching to the 
scattering solutions (details in the Appendix) yields the t- 
matrix. Then the i^ a6 (e, r) irregular wavefunction is cal¬ 
culated similarly by an integration inwards starting at the 
muffin-tin radius. The integrations are performed with a 
predictor-corrector algorithm 26 on logarithmic scale with 
721 radial mesh points in the muffin-tin sphere. To ob¬ 
tain the normal self-consistent potential, V e ff(r), the 
energy integrals are performed by sampling 16 points 
on a semicircle contour in the upper complex energy 
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plane. The calculations are carried out within the atomic 
sphere approximation with an angular momentum cutoff 
of Imax — 2. We use 2450 k points for integration over 
the Brillouin zone to calculate the DOS of bulk Nb. 

In what follows, we calculate the DOS and the Bloch 
spectral function (BSF) which is equivalent to the quasi¬ 
particle spectrum. In all of the following plots the energy 
is measured in units of Rydbergs and k in units of tt/cl. 
The contour plots of the spectral functions are calculated 
in 400 energy points x 265 ^-points. 


IV. RESULTS 


The BSF is defined as A B (e, k) = ^( e — e n (fc)) and 

can be calculated directly from the Green function. In a 
layered system for layer /, this can be expressed as 

A;||) = —Im Tr Gjj(e,r,k ||). (35) 


Since the BSF is equivalent to the quasiparticle spectrum, 
drawing a contour plot of the BSF as a function of energy 
along specified directions of k is a powerful tool to visual¬ 
ize the quasiparticle states. In a layered system this can 
be done for each layer, based on Eq. (35). The spectral 
functions were calculated by adding a small imaginary 
part of 0.0005 Ry to the energy. 


A. DOS of bulk Nb 

To test our procedures, and to show the effect of the 
A pairing potential on a bulk sytem (Nb), we first per¬ 
formed calculations for the case of bulk niobium using 
the values Ajsfb = 0 Ry and A Nb = 0.01 Ry. The DOS 
can be calculated from the BSF D(e) = f A B (e, fc||) dfcj |. 
The particle-hole symmetry implies that the density of 
the hole-like states are just the reflection of the density 
of electron-like states to the Fermi energy. This is indeed 
the case in our calculations, as it can be seen in Fig. [l] for 
the case of A^b = 0 Ry (left panel). If A^b is nonzero, 
a gap appears around the Fermi level and the size of the 
gap equals to the value of A 


B. Normal state band structure of Nb/Au 
heterostructures 

To demonstrate the power of our new theory for an in¬ 
homogeneous system, we apply it to study the system of 
Au overlayers on Nb(100). Foremost, we made calcula¬ 
tions for the normal state, for two reasons. First, to ob¬ 
tain self-consistent potentials and work functions for the 
BdG calculations. Second, it is important to explore the 
features in the normal state electronic structure to later 
understand the quasiparticle spectrum we are planning 
to calculate. Therefore, self-consistent calculations were 
performed for systems containing a semi-infinite Nb, an 



FIG. 1. (Color online) DOS (arbitrary units) of bulk Nb 
(. Ef = 0.713 Ry) in the case of Ajv& = 0 Ry (left panel) 
and A Nb = 0.01 Ry (right panel). The blue line corresponds 
to the density of electron-like states and the red to the density 
of hole-like states. 


additional 6 Nb layers and subsequently 3, 9, 24, and 93 
Au layers. In Fig. [2] we show the contour plot of the BSF 
for a layer that we considered to be in the ” middle” of the 
appropriate sample and for a layer in the bulk Nb (seen 
in Fig. [3] top left panel). It should be mentioned that the 
later is just the projection of the bulk spectral function 
on the (100) plane and it represents the corresponding 
projection of the bulk band structure. The plots are re¬ 
stricted in energy to the range of [-0.05 Ry,0.05 Ry] (later 
we will choose Ajyb to equal this value, and solve the BdG 
equations within this energy range). When Fig. [ 2 ] viewed 
as a sequence, one can immediately recognize the signa¬ 
tures of confinement. Where the DOS in the bulk Nb is 
low, the states in the Au are confined, as they can not 
scatter into the Nb, and on the other side the system is 
limited by vacuum. In regions where the DOS is high 
in the Nb, the states in the Au are smeared out, as here 
the appropriate electrons can scatter more easily into the 
other side of the interface. The confined states in the Au, 
therefore, can be regarded as quantum-well (QW) states. 
The confinement causes a roughly 2tt/L sampling (where 
L is the thickness of the Au sample) of the Au band 
structure. It can be seen from the figure, that as L in¬ 
creases, the QW bands become denser and denser. As L 
approaches infinity, the bulk electronic structure of Au is 
recovered in the middle of the sample. 


For a fixed number of Au layers, one can investigate 
the layer dependence of the electronic structure. This 
is illustrated for the system with 9 Au layers in Fig. [3] 
First, we have to notice that the QW bands do extend 
into the self-consistent Nb layers, as these layers show 
signatures of both of the bulk Nb and the confined states 
of Au. Surprisingly, around the actual interface there is 
a very sharp horizontal band, which can be seen only at 
the interface layers and quickly disappears further away. 
It is not present either in the bulk Nb, nor in the Au 
electronic structure, therefore, it may be regarded as an 
interface state. 
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Contour plot of A(k,e) for layer 8 (Au) 



Contour plot of A(k,e) for layer 9 (Au) 




FIG. 2. (Color online) Contour plot of the BSF (normal state 
band structure) from the ” middle” of the Au layers for differ¬ 
ent thicknesses of the Au: 3 Au layers (top left panel) and 9 
Au layers (top right panel) 24 Au layers (bottom left panel) 
and 93 Au layers (bottom right panel). 


Contour plot of A(k,e) for bulk Nb 
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Contour plot of A(k,e) for layer 6 (Nb) 
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FIG. 3. (Color online) Contour plot of the BSF (normal state 
band structure) for different layers. In the interface region 
there are 6 Nb layers, 9 Au layers, 3 vacuum layers. 


C. The quasiparticle spectrum of Nb/Au 
heterostructures 

We now consider the solution of the BdG equations 
described in the theory section, we model the pair po¬ 
tential in the inhomogeneous Nb(100)/Au(100) system 
by assigning a constant value A Nb = 0-05 Ry to the Nb 
layers and a constant A a u = 0 Ry for the Au layers (and 
the same to the empty sphere and vacuum layers). The 
results of the calculation are shown in Fig. [4j Similarly 
to the normal state, first, we show results for layers in 
the middle of the systems considered. However, we do 
not show any result for the bulk spectral function, be¬ 
cause it is exactly zero in the energy range of the bulk 


superconducting gap. What can be seen immediately is 
that there is a superconducting gap even in the Au layers. 
This gap must have been induced by the vicinity of the 
Nb, because Aa u = 0 Ry. Examining the details of the 
quasiparticle spectrum, especially the one corresponding 
to the sample with 9 and 24 Au layers, reveals that not 
only one, but in fact several gaps are opened. This is 
in strong contrast to bulk superconductors, where the 
quasiparticle states can be obtained from the electronic 
ones by mirroring them to the Fermi energy and opening 
up a gap. Our result modifies this picture so, that the 
proximity of a superconductor in the studied heterostruc¬ 
tures induces the mirroring of the electronic bands, and 
opens up a gap - which is significantly smaller than the 
one in the bulk - at each band crossing. This is valid 
for those band crossings as well that are not directly at 
the Fermi level but within the Aenergy range. In the 
case of the Nb/Au system - due to the the QW states in 
the normal state - the result is a sort of oscillating quasi¬ 
band in Fig. [4j This is a speciality of the Nb/Au system, 
other systems may not look so clean. The induced gap 
is opened between the mirrored branches of the interface 
state as well. However, in contrast to the QW states, the 
interface states do shift quite significantly upwards in en¬ 
ergy, and these states still disappear rather quickly away 
from the interface. It should be noted as well that those 
regions of the spectrum which were more or less smeared 
out in the normal state, now sharpened up. This is the 
consequence of the opening of the superconducting gap 
in the Nb: the states where scattering into was allowed 
before, now disappeared. 


Contour plot of A(k,e) for layer 8 (Au) w. BdG 


Contour plot of A(k,e) for layer 9 (Au) w. BdG 
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Contour plot of A(k,e) for layer 18 (Au) w. BdG Contour plot of A(k,e) for layer 53 (Au) w. BdG 
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FIG. 4. (Color online) Contour plot of the BSF (quasiparticle 
spectrum calculated from BdG equations) from the ’’middle” 
of the Au layers for different thicknesses of the Au: 3 Au 
layers (top left panel) and 9 Au layers (top right panel) 24 
Au layers (bottom left panel) and 93 Au layers (bottom right 
panel). 

In a quasiclassical picture one expect dispersionless 
ABS. However, this is not what the quasiparticle bands 
are showing. What we find is, that the dispersion of the 



























^-dependent ABS can be understood from the features of 
QW states in the normal state as it was described above. 
In conventional superconductivity, the gap is assumed to 
be ^-independent, while in our calculations the obtained 
energy gap strongly depends on the two-dimensional k. 
This is quite surprising, considering the fact that our 
calculations involved only a totally conventional super¬ 
conductivity scenario. This is in an even larger con¬ 
trast to the result of Suvasini et alP who obtained only 
a very week /^-dependence of the gap in bulk Nb. In this 
sense our results show similarities between the physics of 
conventional superconductor - normal metal heterostruc¬ 
tures and unconventional superconductivity. 


Further interesting features of the quasiparticle spec¬ 
trum are revealed if we analyze the spectrum layer by 
layer for a fixed system size (6 Nb and 9 Au layers, see 
Fig.§. As the QW states did overlap with the Nb lay¬ 
ers in the normal state, they still do in superconducting 
state. However, as the quasiparticle states in the Au show 
a much smaller gap than the one in the Nb, these overlap¬ 
ping states lessen the gap in the Nb layers next to the Au 
interface. By performing further calculations, where the 
interfacial Nb layers were more numerous, we found that 
this effect decays quickly, but can be observed up to 15 
layers. In the other side of the interface, in the Au layers 
the induced gap remains constant for each layer. There¬ 
fore, an induced superconductivity may be observed in 
the Au overlayers. This is in accord of the experimental 
observations, where it was found that the whole Nb/Au 
system is superconducting (a common T c has been ob¬ 
tained experimentally in Refs. |10] and HTjl . Cooper pairs 
can be found in the whole system, and the induced gap 
- that appeared in the quasiparticle spectrum in each of 
the Au layers - can be interpreted as a consequence of 
an effective electron-phonon coupling in the Au overlay¬ 
ers caused by the semi-infinite Nb. Quite surprisingly, in 
our calculations we did not find any layer dependence of 
the induced gap. This can be attributed to the fact that 
we did not consider the layer dependence of the pairing 
potential, or by other words, the layer dependence of the 
electron-phonon interaction was neglected. Nevertheless, 
the size of the gap does change with the thickness of the 
system, as it can be seen in Fig. [5] and also summarized 
in Fig. [6] It shows a fast decay, however, it can not be 
fitted well by an exponential function. 


It is also useful to mention that for k = 0 the spec¬ 
trum of the ABS is comparable with the results of one¬ 
dimensional model calculations and the Andreev energy 
levels show the similar 1/L dependence which were also 
obtained in Ref. ns However, we emphasize that this 
property is the consequence of the roughly 2ix/L sampling 
connected to the QW states and can not be regarded as 
an universal feature for every S/N heterostructures. 


Contour plot of A(k,e) for layer 3 (Nb) w. BdG 


Contour plot of A(k,e) for layer 6 (Nb) w. BdG 
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Contour plot of A(k,e) for layer 12 (Au) w. BdG 


Contour plot of A(k,e) for layer 16 (Vac) w. BdG 
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FIG. 5. (Color online) Contour plot of the BSF (quasiparticle 
spectrum calculated from BdG equations) for different layers 
in the system consisting of 6 Nb layers, 9 Au layers and 3 
empty sphere layers. 



FIG. 6. Induced gap on the Au layers as function of the 
thickness of the Au, extracted from results similar to that 
shown in Fig. [4] 


D. Surface states 

Metallic surfaces often exhibit a Shockley-type surface 
state. The energy of such states are located in a rela¬ 
tive band gap of the bulk, normal state band structure 
and usually have a parabolic dispersion, and therefore 
such electrons behave like a nearly free two dimensional 
electron gas. Surface states are easily accessible to spec¬ 
troscopy with photoemission, since they are often located 
near the Fermi energy. Therefore, it is interesting to 
study such surface states once the material becomes su¬ 
perconducting. 

We calculated such Shockley-type surface state from 
the BdG equations in the case of the investigated Nb/Au 
heterostructure along the direction k y = 0. It should be 
emphasized that this surface state is entirely fictitious, 
as this surface is of an Au(100) in the BCC lattice struc¬ 
ture. First, setting the = 0 Ry (see the left panel of 
Fig.0 , the surface state can be observed in the normal 
state electronic structure. While applying a finite A^b 
pairing potential does open up a gap in the Au, just as 
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we discussed earlier for the case of k x = k y direction, but 
no gap opens at the crossing of the surface state bands, 
indicating that it does not couple to the superconductor. 
This effect can be attributed to the fact that obviously 
the surface state is quite localized to the top layers of 
the metal surface and it is mainly isolated from the bulk 
states. Consequently, they do not take part in the An¬ 
dreev scattering process and thus they do not have a gap 
in the spectrum, as it can be seen in Fig. [7| 

As we indicated earlier, an opposite behavior could 
be observed for the interface state, which is localized to 
the Nb/Au interface. The energy of these states shifts 
upwards. This can be explained by the stronger interac¬ 
tion between the superconductor and the normal metal 
resulted in a larger gap than in the QW states. 


Contour plot of A(k,e) for layer 15 (Au) w. BdG 



Contour plot of A(k,e) for layer 15 (Au) w. BdG 
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FIG. 7. (Color online) Contour plot of the BSF in the k y — 0 
direction corresponding the last layer of Au. The Au sample 
consisted 9 layers. The quasiparticle spectrum was calculated 
from BdG equations. Ajv& = 0 Ry is used on the left panel 
and A Nb — 0.05 Ry on the right panel. 



V. SUMMARY 

In this paper we have presented the first material spe¬ 
cific calculations for an s-wave superconductor - normal 
metal heterostructure. Based on first principles BdG 
equations a novel and unique computer code was devel¬ 
oped which allows us to study the nature of the Andreev 
bound states related to the proximity effect in normal 
metal - superconductor heterostructures. 

For the first time we have extended the SKKR method 
for the solution of the KSBdG Eqs. 0 which allows 
one to investigate the quasiparticle spectrum of super¬ 


conducting heterostructures. In order to compare our re¬ 
sults with normal state electronic structure calculations, 
a scalar relativistic generalization of the BdG equations 
within Multiple Scattering Theory was also provided. 
Formally, the generalized Faulkner-Stocks formula, given 
by Eq. (30), is the main result of this paper. 

To illustrate the power of the new method, it was ap¬ 
plied to Nb/Au heterostructures. For simplicity, Au over¬ 
layers of BCC(100) lattice structure on a Nb BCC(100) 
host have been investigated. While such material is not 
likely to exist for larger Au thicknesses, by assuming a 
layer by layer growth, it resulted in an easily understand¬ 
able system with quantum well states. The effect of the 
superconducting host on the quasiparticle spectrum of 
Au overlayers can be more easily identified by these states 
than on a more complex band structure of a real mate¬ 
rial. Calculations for a more realistic geometry will be 
published elsewhere. 

We showed that the QW states (we found to exist in 
the normal state band structure calculations) become 
bound Andreev states due to Andreev scattering. The 
major result of our investigations is that the ABS have 
dispersion, which can be obtained only by developing the 
BdG-SKKR method. We also found that the proximity of 
a superconductor in the studied heterostructures induces 
the mirroring of the electronic bands, and opens up a 
gap at each band crossing, and the gaps are strongly /in¬ 
dependent. We have seen that this induced gap remains 
constant for each layer for a given Au thickness, however, 
the size of the gap decays as function of the Au thick¬ 
ness. For k = 0, the one-dime_nsional model calculations 
of the Andreev energy levels 19 are recovered for those 
heterostructures where the nearly free electron approach 
is applicable. We also investigated the properties of the 
surface state at the Au surface and found that the gap 
does not appear in the energy spectrum of these states, 
probably, because they are localized to the surface and 
consequently do not take part in the Andreev scattering 
process. 
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Appendix: Mathcing the wavefunctions 

The boundary conditions, at the muffin-tin sphere boundary, can be expressed as follows (a = e, h) for the radial 
part of the wavefunction: 


r m tR?(r = r mt ) = Af P“’ (1) (x = x mt ) + Bf P^ 2 \x = x mt ), 


r m t ^ (' rR?(r )) 


4° ^ a ’ (1 U) 


+ b i ^ p r ,(2 U) 


(A.l) 

(A.2) 
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where P^ ,(yl \x) and Pj* )(y2 \x) are the regular solutions of the scalar relativistic BdG Eqs. (19) inside the muffin-tin 


sphere, and r m * is the radius of the muffin-tin sphere (x m t = logr m t). We emphasize that there are two independent 
regular and two independent irregular solutions of Eqs. ( [l9| ). 

The matching conditions can be written in matrix form: 


M a e = b e , M a h = b h , 


where 



(P 

Bf 

h 

(Al\ 

Bj 1 

a e = 

tf 

, Ql = 

tf h 


P) 


Vi h ) 


(A.3) 


(A.4) 


b e = 


/ T‘mtjl(P Tmt ) ^ 

0 

rmtjl(p e rmt ) + r mtP e jl(p e rmt) 

V 0 ) 


b n = 


( 0 \ 

rmtji{p h rmt ) 

0 

\ r mtjl(p hr mt ) + r mtP h jl(p hr mt)J 


(A.5) 


M = 


Pp\x mt ) 

(pmt) 


PP) 


^,(2) 

2) 


( pmt ) 
( pmt ) 


ip 6 Gnt/h + (p e Gnt) 


0 


0 

0 


\ 


d x P l e ’ < ' 1 \x mt ) d x Pp 2 \x mt ) ir mt p e (l + p e r mt d r )h'f(p e r mt ) 

\d x Pi' {1) {x mt ) d x Pp\x mt ) 0 -ir mt p h (l+p h r mt d r )hf(p h r mt )J 


l ^ x ± i 

The regular wavefunctions can be continued inside the muffin-tin sphere as follows 
• for electron-like incoming wave 


r (R?{r) 

{R? e (r) 


ji{p e r ) - i p e t e l e h'l(p e r) 

lp h fheh+ 


A? 


Pp\r) 


PP) 


(r). 


+ Bf 


Pp\r) 


p->( 2 ) 


(r) t 


(A.6) 


(A.7) 


• for hole-like incoming wave 


- (R?(r) 

W(r) 


= r 


— i if'tf'hl ( p e r ) 
ji(p h r ) + ( p h r ) 


A, 


D e >(i) 

3 ^( 1 ) 


(r) 


+ £/ 


3 e,(2) 

A (2) 


(r) 

C), 


(A.8) 


and the irregular solutions as 


jee( r ) q 

0 J/^(r) 


—>► 


/f e (r) If h {r) 
Ip(r) I? h {r) 


(A.9) 


Also, to calculate the Green function, given by Eq. (|30j), the determination of the normalized irregular solution, in¬ 
side the muffin-tin sphere, is indispensable. To satisfy the matching conditions, one needs to use the linear combination 
of the regular solutions, P^ 1 \x), P^^ 2 \x), and the irregular solutions, P^^ 1 \x), Pj*^ 2 \x): 


I? e {x) 

ipw 


= Af 


3 e,(l) 

l 

PP) 


(x) 

( x ). 


+ Bf 


PP 2 \X)\ d e 


Se,(l) 

P->(X) ( 


(x) 


■ Ppix), 


Df 


5 e i(2) 

A. (2) 


( x ) 

(x) t 


(A.10) 


If h {x) 

ifHx) 


= A f 


3 e,(l) 

l 

PP) 


(x) 


b 1 ; 


( 2 ) 
A (2) 


(*) 

(x)j 




D e P) 

l 

PXP 


0*0 

0 * 0 , 


D, 


A,(2) 
A. (2) 


(*) 

0*0, 


(A.ll) 


d e 

VA e / 


= N 


-1 


Xmtjl^P I'mt) 

0 

r m t (1 + P e r m td r ) ji(p e r mt ) 

V 0 


(A.12) 


where 
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( A l\ 

( o \ 

Bi 

= N -1 

rmtjl{p Tmt) 

Cf 


0 

\E) 


ymt (l Y P ^mt^r) jl(.P ^mt) J 


d x Pi' (1 \x mt ) 

\d x p?' {1 \x mt ) 


Pi e ' {2 \xmt) 

Pl' { 2 \xmt) 

d x Pl' {2 \xmt) 

d x Pi' {2 \x mt ) 


Pl’ (1 \xmt) 

P^M 

d x Fj’ W (x mt ) 

d x Pi' {1 \x mt ) 


Pl' {2 \xmt) \ 
P^ 2 \x mt ) 
d x Pf’ ( - 2 \x mt ) 
d x P^ 2 \x mt )J 


(A.13) 


(A.14) 
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